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Abstract 

A multidimesional function y{r) defined by a sample of points {fiji/i} is approx- 
imated by a differentiable function y(r). The problem is solved by using the Gauss- 
Hermite folding method developed in the nuclear shell correction method by Strutinsky. 

1 Introduction 

Our aim is to approximate a sample of points {fi, yi} which represents a measured or hard 
to evaluate data by a differentiable function y(x). We would like to solve this problem using 
the Gauss-Hermite folding method which idea was originally proposed by V.M. Strutinsky 
PP and later-on generalized in Ref. 0. A detailed description of this method may be found 
already in text-books, e.g. in Ref. [Sj. Having the width of the folding function comparable 
with the average distance between Xi points in the i*^ direction one can obtain the folded 
function which goes very close to the data points but increasing its width one can also wash 
out the fine structure stored in the data. Usually the Strutinsky method was used to realize 
the second scope. The parameter of the folding procedure will be determined by requirement 
that the integral in the i^^ direction of the folded function should be the same as the integral 
evaluated with {xi, yi} pairs using the trapezium rule. A corresponding Fortran program for 
the approximation in the ra- dimensional space is listed in Appendix. 

2 General folding formulae in the one-dimesional case 

We consider an ensemble of N points {xi} distributed uniformly in the interval [a, 6]^. To 
each point Xi corresponds a point y^, and we assume there exists a function y{x) such that: 

Vi = y{xi) . (1) 

Let jn{x,x') be a symmetric function of its arguments (i.e. jn{x,x') = jn{x',x)) having 
the following properties: 

+ 00 

J jn{x,x') dx = l (2) 

— oo 

and 

+ 00 

Pfc(x) = / Pk{x')jn{x,x')dx' , (3) 



^ Truly speaking the assumption about the uniform distribution of the points is too strong. It is sufficient 
to assume that the points Xi have to cover the whole interval [a, b] and to be ordered i.e. Xi+i > xi). 
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where k < n are the even natural numbers and Pk is an arbitrary polynomial of the order 
k. In the following the function jn{x,x') will be called the folding function of the n*^ order. 
An example of such a folding function can be a combination of the Gauss function and the 
Hermite polynomials of argument proportional to |x — x'| frequently used in the Strutinsky 
shell correction method [H Ej . More detailed description of the Strutinsky folding function 
will be given in the next section. 

With each discrete point one can associate the function yi{x) defined by: 

+00 

Viix) = / yiS{x' - Xi) jn{x,x')dx' , (4) 



where 6{x) is the Dirac 5-function, so that 

Viix) = yiin{x,Xi) . (5) 

Using Eq. Q it is easy to verify that the integral of the function yi{x) is 

+00 

yi{x) dx = yi . (6) 



Let us construct the function y{x) by summing up all functions yi{x) corresponding to each 
Xi point 

N 

y{x) = ^^iVii^) ■ (7) 

i=l 

The Lebesgue theorem says that the function y{x) is an approximation of y{x) if the weights 
uji are determined from the assumption that the integrals of the unfolded and folded functions 
are (nearly) equal: 



+00 ^ 



/ y{x) dx ^ / y{x) dx = ( 
^ ^ i=i 

a —00 

The Riemann formula for the integral of the function y{x) between bounds a and b reads: 



i=l 



b AT 



/y{x) dx = lim > vixi) A 
i=l 



Xi , (9) 



where Axi is set to: 

Axi = - {xi+i - Xi^i) (10) 

with xq = a and xn+i = b. Comparing Eqs. (jHl) and (jHI) one can see that a reasonable choice 
of the weight is 

Ui = Axi . (11) 
If the number of sample points {xi, yi) large enough than the condition (jH)) will be fuUfiled. 

So, finally the folded function y{x) is given by 

TV 

yi.x) = ^y^Ax^in{x,Xi) . (12) 



1=1 
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3 Gauss-Hermite folding function 



Let the folding function x') be a modified Gauss function 



1 I / / \ 2 I / / 

I I / rr> rr> \ I / />' rr>' 



x') = ^ exp - j /„ — ) , (13) 



where 7 is the parameter and fn{^T~') is the so called correction polynomial of the order n 
determined by the Strutinsky condition In the following we would like to evaluate the 
coefficients of the correction polynomial using some properties of the Hermite polynomials 
which are orthogonal with the weight equal to the Gauss function. 

Let us introduce the variable u = {x — x')/'~f which belongs to the interval (—00, +00). 
The smearing function jn{x,x') and the polynomial Pn{x) (jS)) can be now written as 



-m2 

e 



jn{x, X') = ^ fn{u) , (14) 



Pnix') = Pn{x -1U) = Pn'iu) , (15) 

and 

P„(x) = P„(x + 70) = P„'(0). (16) 
Let us decompose the function Pn{u) into series of the Hermite polynomials Hi{u) 

n 

Pn'iu) = J2^^H,{u) . (17) 

i=l 

Now the condition (jHl) for k = n can be written as 

+ 00 

P„'(0) = -^ j Pn'{u)e~^" Uu)du (18) 

— oo 

and inserting the relation (fTTj) one obtains 



+ 00 
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e-" H,{u) Uu) du - H,{0) \=0. (19) 



oo 



The last equation should be fuUfiled for arbitrary values of 7^ 0) what leads to the following 
set of equations 

/ e-"' H,{u) Uu) du = H,,{0) , (20) 

where i = 0,2, ...,n. From the other side the correction function /„(u) can be also decom- 
posed into series of the Hermite polynomials 

n 

Uu) = Y,CkHu{u) . (21) 

fc=i 

Inserting the above relation to Eq. (j^UI) one obtains 

Hi{0) = V / e-^"H,{u) Hkiu) du . (22) 



Then using the orthogonahty properties of the Hermite polynomials 

"+00 



e-'''Hi{u)Hkiu)du = 2U\Sik , (23) 



one obtains the coefficients of the correction polynomial (PT|) 



C^ = ^Hm (24) 



The values of the Hermite polynomials at zero-point are 



1 for i = , 

i7.(0) = <( 2"(-l)"(2n- 1)!! for i = 2n , (25) 

for i = 2n + 1 , 

so 

r 1 for i = , 

C^ = I (-l)"$(2i^ for z = 2n > , (26) 

[ for 2 = 2ra + 1 . 

The ffist few coefficients Cj and the corresponding Hermite polynomials are: 
Co = 1 Ho = l 



C2 = -\ H^(u) = Au^-2 
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Ci = +^ Hi{u) = I6u^ - 48^2 + 12 

Cg = -gij He{u) = 64^6 - 480m^ + 720^2 - 120 

and the corresponding correction polynomials have the following form 

fo{u) = 1 , 

u{u) = i - f + - , 

Finally the function y{x) approximated by the Gauss-Hermite folding reads: 



(27) 



(28) 



y(a:) = ;^E^.Ax.exp|-^:--:^J j fn ) • (29) 

As a rule the smearing parameter 7 is arbitrary and it can be different at each point Xj. 
But it should be related to the distance Axi between subsequent points if one would like to 
approximate the function stored in the mesh of {xi,yi} points. Similarly one has to choose 
the 7 parameter of the order of the period-length of the fine structure (e.g. shell effects) in 
case when one would like to wash out this structure from the function y{x). 



2 

rr> ly . \ I / 'y nf . 



4 




Figure 1: The sinus function and its approximation with the 2"'^ order Gauss-Hermite 
smoothing. 

4 Multidimensional case 

The extension of the formahsm described in the previous sections is straightforward. Let 
assume that the data points are stored in a m-dimensional array Y[l : Ni,l : A^2, 1 : ^m] 
which corresponds to the ordinates given by the m one-dimensional arrays Xi [1 : Ni] , where 
i = 1,2, ...,m. It means that each element of Y[ii,i2, ....,im] corresponds to the coordinate 
^fc[^fc] with k = 1,2, m. 

The ensemble of the Hermite polynomials Hi{x) forms a complete basis of orthogonal 
functions in which an arbitrary m-dimensional function F{r) = F{xi,X2, ■■■,Xm) can be 
expanded 



oo oo 



F{xi,X2,...,x„,) = ^ ^....^Ci^i2,„i^Hi^{xi)Hi.^{x2) ...Hi^{xm) ■ (30) 



«1=0 42=0 im=0 



The same is true for any polynomial Pn{xi,X2, Xm) of the order n but in this case the 
upper limit of the sums in equation analogous to will be n. It means that the folding 
function in the m-dimensional space is simply the product of the m-th one-dimensional 
foldings performed in each single direction: 

m 

<-^n('^l ) -^li •^2i 2'2i 1 •^m) Jn (-^i) -^i) ■ (31) 

i=l 

The folded function Y{xi,X2, ■■■,Xm) is given by the equation analogous to (fT^: 



(32) 



^ oo oo oo 

Y{Xi,X2, ...,Xm)= E ^^ii E AXi^... E ^^ir.- 

il=0 12=0 im=0 

■ Y[ii,i2, ... ,im]jniXi,Xi[ii]) .. 

The folding function jn{xi,x'j) for the Gauss-Hermite smoothing (fT!^ is 

Ux,,x',) = ^exp\-(^^y\fJ^^'] . (33) 
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The Gaussian width 7j can be different in each coordinate. The correction polynomial /„ 
was already given by Eq. (j2T} and it is tabulated in p8|l for n=0,2,4, and 6. 

5 Some data illustrating the quality of the method 

The second order {n = 2) Gauss-Hermite folding in a four-dimensional space is used. Taking 
into account higher order correction polynomials ()28j) one does not increase significantly the 
quality of the approximation of the function in the middle of the data region but it would 
need a more careful treatment of the problem at the edges. The folding is performed using 
the p mesh points closest to the given point in each direction. The tested function are 
spanned on 21^ points. In addition it is assumed that outside the the data region {xi, i/i} 
the function which should be folded has a constant value (equal to the value of the first or 
the last point in the given direction depending from which side of the data region one takes 
the data for folding). 

The data in Table ^ are listed for some values of the smearing parameter 7 in order to 
see its influence on the accuracy of the approximation. The cosines function in the four- 
dimensional space is chosen as the test function: 

Y{xi, X2, X3, X4) = cos{r) , (34) 

where r = -^/xf + + x| + a;| and the equidistant points Xn{i) G (— 27r, 27r) for i = 1, 2, 21 
and n = 1, 2, 3, 4. The upper and lower limit are Y^ax = 1 and Ymin = —1 respectively. 

Table 1: The approximation of the four- dimensional cos(r) function by the 2"^ order Gauss- 
Hermite folding on basis of p = 5 and p = 7 closest to the given points in function of the 
smearing parameter 7. 



1/7 


^avr 


'^min 


'^max 


^avr 




•^max 


p = 5 


p = 7 


0.98 


0.0081 


-0.0296 


0.0530 


0.0030 


-0.0073 


0.0241 


1.00 


0.0072 


-0.0261 


0.0485 


0.0029 


-0.0074 


0.0242 


1.02 


0.0065 


-0.0233 


0.0452 


0.0030 


-0.0076 


0.0249 


1.04 


0.0060 


-0.0210 


0.0428 


0.0032 


-0.0080 


0.0260 


1.06 


0.0057 


-0.0192 


0.0414 


0.0035 


-0.0086 


0.0276 


1.08 


0.0057 


-0.0179 


0.0409 


0.0040 


-0.0093 


0.0295 


1.10 


0.0059 


-0.0171 


0.0411 


0.0046 


-0.0102 


0.0320 



The root mean square deviation 



^hym]-Y[rm 

i=l 



N 



\ 



2\ 
I 



1/2 



(35) 



as well as the maximal in plus difference (5max) and the minimal in minus one (5min) are 
evaluated for the = 149057 mesh and inter-mesh points (in the middle) excluding the 
points which lie on two outer layers (i.e two first or last rows or columns). Such a choice of 
the test nodes was made in order to eliminate the influence of the border condition on the 
deviation 5»^r- 
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Some other examples of the accuracy of the 2 order Gauss- Hermite approximations are 
hsted in Table |21 The function written in the first column are tabulated at 21 equidistant 
points in the each coordinate in the 4-dimensional space in the interval written in the 2"^^ 
column. The smearing parameter 7 = 0.93 or 7 = 1 is chosen in case of the j9 = 5 or 
p = 7 point basis used when folding, respectively. It is seen in Tables [U and El that in all 

Table 2: A few examples of the approximation accuracy in the four-dimensional space. 



Function Y 


Range 


Y ■ 

-' mm 


Y 

-' max 


P 


•^avr 


•^min 


•^max 


sin(r)/r 


-27r : 27r 


-0.2172 


1 


5 


0.0011 


-0.0029 


0.0128 


7 


0.0005 


-0.0012 


0.0059 


^2 1 ^2 1 ^2 1 /J-.2 


-2 : 2 





8 


5 


0.0053 


-0.0218 


0.0102 


7 


0.0013 


-0.0017 


0.0014 


(Xi ■X2- X3- X4,y 


-2 : 2 





256 


5 


0.0076 


-0.2491 


0.1161 


5 


0.0017 


-0.0191 


0.0102 


Xl ■ X2 ■ X3 ■ Xi 


-2 : 2 


-16 


16 


5 


0.0014 


-0.0180 


0.0180 


7 


0.0001 


-0.0017 


0.0017 



considered cases the root mean square deviation {5 avr) is of the order 10 ^ or less of the 
maximal difference Ymax — Ymin between the data points. 

6 Summary and conclusions 

A new method of the smooth approximation of a function defined on a sample of points 
in a multidimensional space is proposed. The folding of the discrete data points using the 
Gauss-Hermite method of Strutinsky is performed. Depending on the width of the Gauss 
function the folded function can be very close to the approximated data or can give its 
average behavior only. 

The folded function and all its derivatives are continuous. This significantly increases the 
range of applicability of the method. Our approximation of the discrete data can be used e.g. 
when solving transport equations or other type of equations of motion in a multidimensional 
space, what is a frequent problem in economy, meteorology and environment protection 
problems as well as in molecular or nuclear dynamics. 

The proposed approximation of the data can be also used in the computer graphic art. 
It can wash out the fine structure from a photography keeping unchanged its average back- 
ground. One can also think about the use of the new folding method when one evaluates 
the cross-sections of a multidimensional data which one has e.g. in the X-ray tomography. 
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7 Appendix 

The source of the fortran program for the 2"*^ order Gauss- Her mite approximation: 



Q 

C The 2nd order Gauss-Hermite folding (a la V.M. Strutinsky) of a function 
C defined on a sample of equidistant points in the n-dimensional space 

Q 

subroutine ghost (nna, npts, xdn, dx, y, yref , gamma, x, fun, dfun) 
parameter (ndim=4 , id=2 , nid=2 *id+l ) 

dimension y (npts) , xdn (ndim) , dx (ndim) , x (ndim) , dfun (ndim) , nna (ndim) 

C 

C The data points which should be approximated by the function fun (x) have 

C to be stored in the main program as the n-dimensional array: 

C dimension y (nna (1) , nna (2) , . . . . , (nna (ndim) ) , 

C where nna(i) is the number of points related to the x_i coordinate. 

C The total number of points in y is npts=nna (1) *nna (2) * *nna (ndim) . 

C The equidistant grid beginning at xdn(i) with step dx(i) is assumed 

C for each coordinate . The value and gradient of the approximated function 

C in the the point x (ndim) are stored in fun and dfun respectively. 

C REMARK: In order to increase the accuracy of the approximation one 

C preforms the folding of the differences y(i)-yref, where yref is the 

C function value around which the approximation should be the best. 

C A reasonable choice of yref is the average of the input points i.e. 

C yref=(sum y (i) ) /npts . The folding is performed using 2*id+l points 

C closest to point x in each direction and gamma is the smearing width. 

C (C) Copr. 2004 Krzysztof Pomorski, email: pomorski@kft.umcs.lublin.pl 

C 

dimension f (ndim, nid) , df (ndim, nid) , fdf (ndim) , ni (ndim) , nnn (ndim) 

gami=l . / gamma 

f un=0 . 

nnn (1) =1 

do 2 i=l,ndim 

dfun (i) =0 . 

if(i.gt.l) nnn(i)=nnn(i-l) *nna (i-1 ) 

xx= (x (i) -xdn (i) ) /dx (i) +1 . 

ni (i) =int (xx+0 .5) 

f norm=0 . 

do 1 j=l,nid 

t=gami* (xx- (ni (i) +j-id-l) ) 

gauss=exp(-t**2) 

f (i, j) =gami* gauss* (1 . 5-t**2) 

f norm=f norm+f (i, j ) 

1 df (i, j) =gami**2*gauss* (2 . *t**3-5 . *t) /dx (i) 
do 2 j=l,nid 

f (i , j ) =f ( i , j ) /f norm 

2 df (i , j ) =df ( i, j ) /f norm 
nbox=nid* *ndim 

do 5 k=l,nbox 
1 = 1 
f f=l . 

do 3 m=l,ndim 

3 fdf(m)=l. 

icur=k-l 
nn=nid** (ndim-1) 
do 5 i=ndim, 1,-1 
j=icur/nn+l 
ff=ff*f (i, j) 
do 4 m=l,ndim 

if(m.eq.i) f df (m) =f df (m) *df ( i , j ) 

4 if(m.ne.i) f df (m) =f df (m) *f (i, j ) 

l=l+nnn (i) * (min (nna (i) ,max ( 1, ni (i) + j-id-1 ) ) -1 ) 
icur=icur- ( j-1 ) *nn 

5 nn=nn/nid 

fun=fun+ (y (1) -yref) *f f 
do 6 m=l,ndim 

6 dfun (m)=dfun (m) + (y (1) -yref ) *fdf (m) 
f un=f un+yref 

return 
end 
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